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Abstract 

This  report  presents  the  results  of  the  investigation  of  two  important  subjects  associated  with  the 
design  and  certification  in  unitized  airframe  components.  The  first  subject  is  related  to  the  effects 
of  residual  stresses  on  the  structural  stability  of  thin  unitized  components  machined  from 
aluminum  plates,  in  particular  7050-T74  and  7050-T7451  plates.  The  findings  indicate  that 
residual  stresses  introduced  in  a  plate  during  the  rolling  operation  (bulk  stresses)  and  residual 
stresses  introduced  into  a  part  machined  from  the  plate  during  high  speed  machining  should  be 
included  as  a  modeling  consideration  when  designing  thin  unitized  components.  The  second 
subject  is  related  to  the  computation  of  strain  energy  release  rate  in  damaged  laminate  composite 
materials.  Typical  failure  in  the  presence  of  an  initial  defect,  such  as  delamination,  appears  under 
a  mixed  mode  loading,  therefore  it  is  essential  to  have  an  efficient  algorithm  for  the  computation 
of  the  strain  energy  release  associated  with  each  loading  mode  for  the  construction  a  mixed 
mode  failure  criterion  for  the  determination  of  residual  strength  of  unitized  components  made  of 
composite  materials.  The  Virtual  Crack  Closure  Technique  (VCCT)  was  considered  during  the 
Phase  I  project.  It  was  found  that  typical  numerical  implementations  of  the  VCCT  utilizing  the  h- 
version  of  the  finite  element  method  (FEM)  are  unreliable  because  the  results  are  mesh- 
dependent.  A  modification  of  the  method  was  investigated,  involving  a  combination  of  numerical 
and  analytical  computations,  which  is  well  suited  for  its  implementation  with  the  p-version  of  the 
Finite  Element  Method. 

Introduction 

Increasing  emphasis  on  affordability  of  military  systems  has  led  to  a  number  of  advances  in 
airframe  design  and  production.  Unitized  airframe  structural  components  are  replacing  sheet 
metal  built-up  components  to  reduce  part  count  and  assembly  cycle  times  and  costs.  The 
development  of  high  speed  machining  (HSM)  techniques  has  made  it  possible  to  fabricate  thin 
lightweight  structures  that  provide  improved  performance  at  lower  costs.  There  is  a  need  for 
advanced  numerical  methods  for  the  solution  of  problems  associated  with  the  manufacturing  and 
certification  of  unitized  airframe  components.  The  goals  are  to  support  (a)  damage  tolerant 
designs,  (b)  establishment  of  criteria  for  inspection  intervals,  (c)  planning  the  fabrication 
processes  so  that  the  incidence  of  re-working  and  scrapping  of  partially  or  fully  manufactured 
parts  is  substantially  reduced,  and  (d)  evaluation  of  options  for  reworking  and  clear  criteria  for 
decisions  whether  to  certify  or  reject  an  as-manufactured  part.  This  is  expected  to  result  in 
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substantial  affordability  improvements  for  aircraft  and  spacecraft  structures  as  well  as  increased 
reliability. 

Most  aluminum  and  many  titanium  aircraft  components  are  fabricated  by  milling  out  90  to  95 
percent  of  the  material  from  plate  stock.  Typical  parts  consist  of  large  numbers  of  open  pockets 
with  thin  bottoms  (webs),  enclosed  by  thin  walls  (ribs)  perpendicular  to  the  webs.  The  webs  and 
ribs  are  connected  by  filleted  transition  regions.  The  problem  class  is  sufficiently  large,  important 
and  complex  to  warrant  development  of  specialized  software  for  analysis  and  design. 

It  is  necessary  to  have  the  capability  to  verify  and  validate  the  mathematical  models  that  serve  as 
a  basis  for  engineering  decisions.  Parametric  modeling  capabilities  based  on  the  hierarchic 
concept  of  models,  that  allow  dimensional  reduction  where  appropriate,  and  account  for  non- 
linearities  when  necessary,  provide  the  mathematical  and  technological  basis  for  reliable 
numerical  simulation  procedures. 

The  problem  of  buckling  (oil  canning)  was  observed  in  machining  experiments  conducted  at  the 
Machining  Development  Laboratory  of  the  Boeing  Phantom  Works,  St.  Louis,  Missouri.  It  was 
noted  that  thin-walled  test  articled  buckled  upon  removal  from  the  fixtures.  The  results  of 
preliminary  investigations  suggest  that  the  machining-induced  residual  stresses  are  the  primary 
cause  for  buckling. 

The  Phase  I  effort  focused  on  the  conceptual  basis  of  an  algorithmic  structure  designed  to  meet 
the  objectives  indicated  above  and  development  of  an  implementation  plan.  ESRD  consulted  with 
AFRL  Materials  and  Manufacturing  Directorate  and  aerospace  OEMs  (The  Boeing  Company  in 
St.  Louis,  Missouri  and  Lockheed  Martin  Aeronautics  in  Marietta,  Georgia)  in  the  formulation  and 
prioritization  of  the  technical  objectives.  During  these  consultations  it  became  apparent  that  a 
useful  tool  for  the  analysis  and  design  of  unitized  structural  component  should  also  include  a 
capability  to  determine  the  residual  strength  of  damaged  structures  made  of  composite  materials. 
The  strain  energy  release  rate  was  identified  as  the  key  parameter  that  characterizes  the  residual 
strength  of  unitized  structures  made  of  laminate  composites. 

Since  composite  materials  generally  exhibit  a  coupling  between  load  and  deformation  modes 
(i.e.,  a  symmetric  loading  may  lead  to  a  deformation  that  is  not  purely  symmetric)  due  to  lack  of 
symmetry  in  material  properties.  Mode  I,  Mode  II  and  Mode  III  loading  are  expected  to  be  present 
in  most  cases  when  fracture  mechanics  parameters  for  composite  materials  have  to  be 
determined.  The  most  promising  characterizing  parameters  for  the  onset  of  propagation  in 
composites  is  the  energy  release  rate  associated  with  each  mode  of  deformation.  During  the 
Phase  I  project  the  virtual  crack  closure  technique  (VCCT)  within  the  framework  of  the  p-version 
of  the  finite  element  method  was  investigated  as  the  algorithm  for  the  computation  of  the  energy 
release  rate  for  Mode  I  and  Mode  II  in  composite  materials  for  2-dimensional  problems. 
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The  following  tasks  were  performed  during  the  Phase  I  investigation: 

1 .  Preliminary  investigation  of  the  effects  of  residual  stresses  caused  by  (a)  the 
manufacturing  process  of  7050-T74  and  7050-T7451  aluminum  plate  stock  and  (b) 
mechanical  milling  on  the  structural  stability  of  unitized  metallic  airframe  components. 
Machining-induced  stresses  result  from  the  interaction  between  the  work  piece  and  the 
cutting  tool.  These  stresses  are  estimated  to  be  of  the  order  of  150  MPa,  within  a 
boundary  layer  of  approximately  0.3  mm.  The  stresses  decay  rapidly  outside  of  the 
boundary  layer.  This  poses  a  challenging  computational  problem,  because  elements  near 
the  boundary  need  to  have  very  large  aspect  ratios.  During  the  Phase  I  investigation  the 
available  residual  stress  profiles  were  utilized  to  determine  the  sensitivity  of  the  classical 
buckling  strength  to  the  magnitude  and  distribution  of  residual  stresses. 

2.  An  algorithm  for  the  determination  of  the  buckling  strength  of  unitized  metallic  airframe 
components  that  accounts  for  the  effects  of  residual  stresses  was  developed.  In  view  of 
the  fact  that  the  loading  data,  material  properties  and,  to  a  lesser  extent,  the  geometric 
properties  are  stochastic,  the  algorithm  needs  to  be  enhanced  to  allow  rapid  evaluation  of 
alternative  designs,  support  parametric  optimization  procedures  and  Monte  Carlo 
simulations. 

3.  Preliminary  investigation  of  the  proper  algorithm  for  the  computation  of  energy  release 
rates  for  defects,  such  as  delamination  in  unitized  structures  made  of  composite 
materials. 

ESRD’s  software  product  StressCheck  was  utilized  during  the  Phase  I  investigation. 

Residual  stresses 

Effect  of  residual  stresses  on  pre-load  buckling 

Consider  the  box  like  test  specimen  machined  from  an  aluminum  plate^  (Figure  1).  Prior  to 
removing  the  specimen  from  the  machining  fixture  no  signs  of  buckling  were  present,  but  after 
removing  the  specimen  from  the  fixture  it  buckled  in  a  typical  ‘oil  can’  mode. 


^  Provided  by  The  Boeing  Company,  Phantom  Works,  St.  Louis  MO. 
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Figure  1:  Test  specimen  on  fixture  frame. 

Since  buckling  is  observed  after  removing  the  specimen  from  the  fixture,  residual  stresses  alone 
can  be  considered  to  be  the  source  of  the  instability  (therefore  the  use  of  the  term  pre-load 
buckling).  Residual  stresses  can  be  classified  into  two  types,  the  first  are  the  residual  stresses 
caused  by  the  metal  forming  process,  called  bulk  or  material  stresses;  the  second  are  the 
residual  stresses  introduced  during  the  machining  process  (i.e.,  machining-induced  residual 
stresses).  To  determine  whether  the  bulk  stresses,  the  machining  induced  stresses  or  a 
combination  of  both  is  causing  the  specimen  to  buckle,  two  sets  of  numerical  studies  were 
performed  using  a  finite  element  representation  of  similar  dimensions  as  the  test  specimen 
(Figure  2).  The  effects  of  one  type  of  residual  stresses  only  were  considered  for  each  set. 


Figure  2:  Test  specimen  and  finite  eiement  mesh. 

For  the  numerical  simulation,  the  fillets  between  the  web  and  the  flanges  were  omitted  from  the 
model.  This  simplification  is  acceptable  for  buckling  analysis  since  the  presence  of  these 
geometric  details  do  not  significantly  affect  the  results. 

Effect  of  bulk  residual  stresses  on  pre-load  buckling 

The  model  shown  in  Figure  2  was  loaded  with  bulk  residual  stresses  typical  of  those  in  76  mm 
thick,  7050-T74  (heat  treated  only)  and  7050-T7451  (heat  treated  and  stretched)  aluminum  plates 
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shown  in  Figure  3  ([1],  [2]).  An  eigen-value  buckling  analysis  was  performed  in  order  to  obtain  the 
critical  buckling  load  and  mode  shape,  for  different  locations  of  the  test  specimen  through  the 
thickness  of  the  plate  as  indicated  schematically  in  Figure  4.  (The  parameter  ex  is  a  measure  of 
the  eccentricity  of  the  box  with  respect  to  the  mid-plane  of  the  plate). 


Depth  (2x3/h) 

a) 

Figure  3:  Typical  Residual  stresses  in  a  76  mm 


b) 

aluminum  plate:  a)  7050-T74  b)  7050-T7451. 


Figure  4:  Location  of  the  box  with  respect  to  the  plate  where  it  is  machined  from. 


The  study  was  restricted  so  that  the  box  was  contained  within  the  interval  -0.  8  <  2X3//?  <  0.8 

because  there  is  no  information  of  residual  stresses  outside  this  interval  (indicated  by  dots  in 
Figure  3)  and  extrapolation  would  not  be  reliable.  The  coordinate  X3  is  in  the  thickness  direction, 
and  h=7Q  mm. 

The  computed  load  factors  and  typical  mode  shapes  are  shown  in  Table  1  and  Figure  5 
respectively.  The  load  factor  is  the  multiplying  factor  that  must  be  applied  to  the  residual  stresses 
to  produce  buckling.  If  a  load  factor  is  greater  than  1.0,  then  the  specimen  will  not  buckle  under 
the  given  residual  stress  distribution. 


Table  1:  Buckling  Load  Factor  for  boxes  machined  from  7050-T74  and  7050-T7451  aluminum  plates. 


Location 

ex 

[mm] 

Load  Factor 
7050-T74 

Load  Factor 
7050-T7451 

1 

30.32 

3.39 

6.12 

2 

21.51 

4.29 

7.28 

3 

12.7 

15.78 

11.78 

4 

3.89 

6.53 

14.64 

5 

-4.92 

3.48 

5.59 
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Figure  5:  Typical  mode  shapes  (1^^  mode)  for  pre-load  buckling  under  bulk  stresses.  Far  left  for  all 
values  of  ex  for  7050-T74.  Center  for  ex=30.32,  and  right  for  ex=3.89  for  7050-T7451. 

The  residual  stresses  were  incorporated  into  the  model  shown  in  Figure  2  by  specifying  formula- 
based  coefficients  of  thermal  expansion  as  described  in  [1].  The  dimensions  of  the  model  were 
selected  such  that  they  are  representative  of  the  specimens  shown  in  Figures  1  and  2 
(76.2x101.62x25.4  mm^,  with  3.81  mm  thick  walls). 

The  investigation  showed  that  the  bulk  residual  stresses  are  insufficient  (and  very  far)  from 
causing  the  specimen  to  buckle  and  that,  in  some  cases,  the  associated  buckling  mode  was  very 
different  from  the  observed  ‘oil  can’  mode  in  the  test  specimens.  It  was  therefore  concluded  that 
the  bulk  stresses  (from  7050-T74  or  7050-T7451  aluminum  plates)  cannot  cause  the  buckling 
observed  in  the  test  specimens.  It  follows  that  machining-induced  residual  stresses  are  the 
probable  cause  for  buckling  of  thin  parts  machined  from  metal  plates.  This  is  discussed  in  the 
following  section. 

Effect  of  machining-induced  residual  stresses  on  pre-load  buckling 

In  the  study  described  in  the  previous  section  it  was  determined  that  the  bulk  stresses  present  in 
7050-T74  or  7050-T7451  aluminum  plates  cannot  be  the  sole  cause  of  buckling  after  machining. 
In  this  section,  results  are  presented  for  the  model  loaded  with  typical  machining-induced  stress 
profiles,  neglecting  bulk  stresses. 

The  machining-induced  residual  stress  profiles  used  in  this  study  were  obtained  from  experiments 
performed  on  specimens  machined  from  7050-T7451  aluminum  plates  with  similar  tools  (end  mill) 
as  those  used  for  analyzing  the  test  specimens  shown  in  Figures  1  and  2.  Figure  6  shows  a 
schematic  of  the  test  specimens  used  for  computing  the  machining-induced  stresses  (see  [2],  [3], 
[4]  for  details).  The  set  of  specimens  consisted  of  two  groups,  one  corresponding  to  specimens 
machined  with  the  side  of  the  tool  (rib  specimens),  and  the  other  corresponding  to  specimens 
machined  with  the  head  of  the  tool  (web  specimens),  as  shown  in  Figure  7.  This  distinction  was 
made  not  only  because  residual  stresses  induced  by  machining  the  material  with  the  head  or  side 
of  the  tool  were  expected  to  be  different,  but  also  because  most  thin  unitized  components  can  be 
considered  as  a  collection  of  ribs  and  webs. 
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X3 


Figure  6:  Test  specimen  description  used  in  the  computation  of  machining-induced  residuai 

stresses. 


a)  b) 

Figure  7:  Test  specimens:  (a)  Web  specimen;  (b)  Rib  specimen. 

Since  experimental  information  of  the  residual  stress  profiles  around  corners  of  fillets  is 
unavailable  these  are  not  considered  in  this  investigation.  It  is  expected  however,  given  the 
rigidity  of  these  components,  that  this  simplifying  assumption  will  not  significantly  affect  the 
results.  Figures  8  and  9  show  the  residual  stress  profiles  used  in  this  study. 


Figure  8:  Residual  stress  profiles  for  a  thin  rib  (ox  and  Oy  respectively) 
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Figure  9:  Residual  stress  profiles  for  a  thin  web  (ox  and  Oy  respectively) 


Machining-induced  stresses  resulting  from  the  interaction  between  the  work  piece  and  the  cutting 
tool  are  of  the  order  of  150  MPa,  within  a  boundary  layer  of  approximately  0.3  mm.  The  stresses 
decay  rapidly  outside  of  the  boundary  layer.  This  poses  a  challenging  computational  problem, 
because  elements  near  the  boundary  need  to  have  very  large  aspect  ratios.  This  difficulty  is 
easily  overcome  by  the  p-version  of  the  finite  element  method  [8].  In  order  to  provide  an  accurate 
approximation,  the  finite  element  mesh  was  designed  with  four  elements  through  the  thickness  of 
the  specimen  as  shown  in  Figure  10.  The  quality  of  the  solution  was  ascertained  by  running  8 
linear  solutions  from  p=1  to  8,  and  checking  the  convergence  of  the  linear  solution  before 
performing  the  buckling  analysis. 


Figure  10:  Mesh  layout  and  detail. 

The  first  buckling  mode  is  shown  in  Figure  11a,  where  the  computed  load  factor  (Ld.  Fct.)  was 
2.34e-01.  It  follows  that  a  specimen  machined  with  the  same  tool  that  produced  the  residual 
stress  profiles  used  in  this  study  will  buckle  once  removed  from  the  fixture.  The  second  buckling 
mode,  shown  in  Figure  11b,  is  the  typical  ‘oil  can’  mode.  The  load  factor  was  computed  as  1.433, 
which  is  very  close  to  1 .  Given  the  fact  that  the  machining-induced  stresses  are  rough  estimates 
only,  it  is  probable  that  this  mode  will  occur  also. 
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Ld.  Fct.=  2.340e-001 


Figure  11a:  First  buckiing  mode  and  ioad  factor. 


Figure  11b:  Second  buckiing  mode  and  ioad  factor. 

At  this  time  no  reliable  methods  exist  for  the  determination  and  mathematical  representation  of 
machining-induced  stresses.  Pre-load  buckling  in  parts  machined  from  aluminum  metal  plates 
may  be  prevented  by  removal  of  the  layer  of  material  affected  by  the  machining  process  (denoted 
as  plastic  layer)  using  a  process  such  as  chemical  milling  (which  is  supposed  to  introduce  no 
residual  stresses)  as  described  in  [2].  By  removing  the  plastic  layer,  the  specimen  will  be  affected 
mainly  by  bulk  stresses,  for  which  it  is  expected  that  pre-load  buckling  will  be  prevented. 

Effect  of  bulk  residual  stress  on  strength  buckling 

In  most  analysis  performed  using  the  theory  of  elasticity,  initial  stresses  are  neglected,  the  main 
reason  being  the  lack  of  available  information  about  their  distribution.  This  is  a  reasonable 
assumption  when  the  part  to  be  analyzed  has  thick  walls.  In  the  case  of  thin  unitized  components 
the  effects  of  initial  stresses  are  more  relevant  and  should  be  considered  in  formulating 
mathematical  modes.  In  the  studies  described  in  the  previous  sections  it  was  concluded  that 
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machining-induced  stress  can  lead  to  structural  instability  in  the  absence  of  external  loads,  where 
the  bulk  or  material  stresses  (present  in  7050-T74  and  7050-T7451  aluminum  plates)  are 
insufficient  to  cause  pre-load  buckling.  It  is  expected  however,  that  the  presence  of  material 
stresses  can  affect  structural  stability  to  a  sufficient  degree  to  cause  a  significant  reduction  in  the 
load  carrying  capacity  of  a  structure. 

In  this  section  we  investigate  the  effect  of  bulk  stresses  on  the  buckling  strength  of  a  simple  test 
specimen  under  a  constant  shearing  load,  as  shown  in  Figure  12. 


Figure  12:  Box  specimen  under  pure  shear. 


The  load  required  to  cause  the  specimen  to  buckle  was  computed  for  three  cases: 

•  Reference  case:  Box  loaded  by  shear  with  no  residual  stresses. 

•  Box  loaded  by  shear  and  including  the  initial  bulk  stresses  of  7050-T74  aluminum  plates 
shown  in  Figure  3a. 

•  Box  loaded  by  shear  and  including  the  initial  bulk  stresses  of  7050-T7451  aluminum 
plates  shown  in  Figure  3b. 

The  computation  of  the  buckling  load  was  performed  by  an  Iterative  procedure  in  which  the 
residual  stresses  were  kept  constant  and  the  shear  load  was  increased  until  the  buckling  load 
factor  was  one.  The  effect  of  residual  stress  in  the  structural  stability  of  the  box  is  represented  by 
the  relative  difference  between  the  reference  solution  and  the  solution  including  initial  stresses.  A 
positive  sign  means  that  the  residual  stresses  are  increasing  the  load  required  to  produce 
structural  instability  (beneficial  effect  of  residual  stresses).  On  the  other  hand,  a  negative  sign 
means  that  the  residual  stresses  are  decreasing  the  load  required  to  produce  structural  instability 
(detrimental  effect  of  residual  stresses). 

Table  2  shows  the  computed  values  for  the  test  specimen.  It  is  observed  that  the  residual 
stresses  affect  mainly  the  buckling  load  whereas  the  buckling  mode  (shown  in  Figure  13)  remains 
unchanged.  The  residual  stresses  present  in  7050-T74  aluminum  plate,  at  a  location  of  ex=30.32 
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mm,  increases  the  stability  of  the  specimen  by  38%,  while  the  residual  stresses  in  7050-T7451  at 
the  same  location  reduces  the  buckling  strength  by  11%. 

Table  2:  Shear  load  required  to  buckle  the  specimen,  with  and  without  bulk  residual  stresses  (RS). 


Buckling  Shear  Load  (MPa) 

ex=30.32 

RS  +  Load 

Load  alone 

Stability  variation 

7050-T74 

0.2272 

0.1645 

38.12% 

7050-T7451 

0.1463 

0.1645 

-11.06% 

Figure  13:  Buckling  mode  for  box  specimen  under  shear  load  (with  and  without  bulk  residual 

stresses). 

This  investigation  indicates  that  the  stability  of  loaded  thin-walled  unitized  structural  components 
depends  not  only  on  the  magnitude  of  the  residual  stresses  but  also  on  the  type  of  loading.  If  the 
compressive  stresses  resulting  from  external  loads  combine  with  the  compressive  residual 
stresses,  then  the  effect  will  be  to  reduce  buckling  strength. 

It  is  well  known  that  substantial  differences  exist  between  predictions  of  buckling  strength  based 
on  mathematical  models  and  the  results  of  physical  experiments.  These  differences  are  usually 
attributed  to  geometric  imperfections  not  incorporated  in  the  model.  The  results  of  our 
investigation  indicate  that  the  effects  of  residual  stresses  are  significant  also. 

Delamination  of  composite  materials 

The  strain  energy  release  rate  in  damaged  laminate  composite  materials  has  been  identified  as 
the  parameter  that  characterizes  the  residual  strength  of  unitized  structures  made  of  laminate 
composites  [5],  Typical  failure  in  the  presence  of  an  initial  defect,  such  as  delamination,  appears 
under  a  mixed  mode  loading.  Therefore  is  essential  to  have  an  efficient  algorithm  for  the 
computation  of  the  strain  energy  release  rates  (G,,  /=1 ,  2,  3)  associated  with  each  loading  mode 
(opening  or  tension  mode  /=1,  sliding  shear  mode  /=2,  and  scissoring  shear  mode  /=3)  for  the 
construction  a  mixed  mode  failure  criterion. 
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The  Virtual  Crack  Closure  Technique  (VCCT)  is  the  most  commonly  used  method  for  the 
computation  of  G,.  The  VCCT  is  based  on  the  relationship  between  the  energy  release  rate  (G) 
and  the  work  required  to  return  a  crack  to  its  original  length.  This  involves  the  computation  of  the 
stresses  in  the  vicinity  of  the  crack  tip  along  a  virtual  crack  extension  of  size  Aa.  Since  the 
stresses  at  the  crack  tip  are  infinity,  it  is  clear  that  in  a  numerical  approximation  this  value  cannot 
be  achieved  (even  when  the  numerical  error  measured  in  energy  norm  can  be  arbitrarily  reduced 
by  using  geometric  mesh  refinement).  Furthermore,  the  singular  behavior  of  the  stresses  at  the 
crack  tip  will  affect  the  solution  near  the  crack  tip.  Therefore  the  accuracy  of  the  computed  value 
of  G  will  deteriorate  when  the  crack  extension  Aa  goes  to  zero.  It  follows  that  the  quality  of  the 
solution  near  the  crack  tip  will  affect  the  computed  strain  energy  release;  therefore  it  should  not 
be  included  in  the  computation  of  the  energy  release  rate  (G). 

It  was  found  during  our  investigation  that  typical  numerical  implementations  of  the  VCCT  utilizing 
the  h-version  of  the  finite  element  method  (FEM)  are  inaccurate  because  the  results  are  mesh- 
dependent.  Typical  applications  used  in  the  h-version  of  the  finite  element  method  utilize  an 
integral  form  of  the  VCCT  equations  which  cancel  some  of  the  effects  of  the  singular  behavior  of 
the  solution  at  the  crack  tip.  However,  the  singularity  at  the  crack  tip  will  affect  the  numerical 
solution  in  the  vicinity  of  the  crack  tip  where  the  computation  of  the  energy  release  rate  becomes 
mesh  dependent. 

A  modified  algorithm  using  the  p-version  was  developed  to  overcome  this  limitation.  This 
algorithm  involves  the  numerical  computation  of  the  stresses  close  to  the  crack  tip,  in  a  region 
where  they  can  be  computed  with  sufficient  accuracy  (guaranteed  by  using  a  geometric  mesh 
together  with  a  uniform  p-extension  in  order  to  ascertain  the  convergence  of  the  solution).  This 
information  is  used  for  computing  the  energy  release  rate  for  a  reduced  crack  path  [e,  Aa],  where 
e  is  an  arbitrarily  small  distance  measured  from  the  crack  tip.  Then,  since  the  behavior  of  the 
exact  solution  in  the  neighborhood  of  the  crack  tip  is  known,  the  computed  solution  is  used  for 
extrapolating  the  stress  field  towards  the  crack  tip  on  the  interval  [0,  e]  using  a  closed-form 
equation,  that  is  used  subsequently  for  obtaining  a  correction  for  the  computed  energy  release 
rate.  In  this  Phase  I  project  we  addressed  the  computation  of  G|  and  Gn  for  two  dimensional 
cases  using  the  p-version  of  the  FEM. 

The  Virtual  Crack  Closure  Technique 
Formulation 

The  Virtual  Crack  Closure  Technique  (VCCT)  is  based  on  the  relationship  between  the  energy 
release  rate  (G)  and  the  work  required  to  return  a  crack  extension  to  its  original  length  (c.f.  [5]  and 
the  references  therein).  A  brief  description  is  included  in  the  following. 
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Consider  the  state  of  stresses  at  the  crack  tip  in  an  isotropic  material  for  a  symmetric  (Mode  I) 
loading.  For  the  coordinate  system  centered  on  the  crack  tip,  as  shown  in  Fig.  1,  and  neglecting 
terms  of  higher  order,  the  Oy  stress  distribution  along  y=0  is  given  by: 

cr  x>0  (1) 

yllm: 

The  displacement  of  the  (top)  crack  face  is: 


u  = 

y 


(l  +  v)(;t-  +  l)  Kj 


x<0 


E  4^ 

Where  E  is  the  modulus  of  elasticity,  v  is  the  Poisson’s  ratio  and  k  is  defined  as: 


(2) 


-  for  plane  stress 

1  +  v 

3  -  4v  for  plane  strain. 


Assume  now  that  the  crack  length  increases  by  a  small  amount  Aa,  as  shown  in  Figure  15.  In  this 
case  the  displacement  will  be: 


(l  +  p)(;t-  +  l)  K, 
E  4^ 


0  <  X <  Aa  ,  x'- x-Aa 


Figure  15:  Crack  tip  coordinate  system.  Notation. 


The  work  required  to  return  the  crack  to  its  original  length,  that  is,  to  close  the  length  increment 
Aa,  is  given  by: 


Aa  A 


—  a„u„dx  - 

2  y  y 


(l  +  v)(/t-  +  l)  K 


(3) 


AW  = 


(l  +  vX/^  +  l) 

4E 


KjAa 


Hence: 


lAanll 
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In  order  to  restore  the  crack  to  its  initial  length,  energy  equal  to  Al/l/  had  to  be  imparted  to  the 
elastic  body.  This  is  the  energy  expended  in  crack  growth,  called  Griffith’s  surface  energy.  The 
potential  energy  had  to  decrease  by  the  same  amount  when  the  crack  increment  occurred. 
Hence: 

^^0  Aa  oa  4E 

Therefore  the  energy  release  G  can  be  computed  as: 


AW 
-  hm - 


Aa 

When  the  loading  is  purely  anti-symmetric  (Mode  II  loading),  then  the  relationship 
energy  release  rate  and  the  work  required  to  return  the  crack  to  its  original  length  is 
the  symmetric  case.  Instead  of  equations  (1)  and  (2)  we  have 


between  the 
analogous  to 


x>0 


The  displacement  of  the  (top)  crack  face  is 


x<0 


(4) 


(5) 


In  view  of  the  fact  that  the  solutions  corresponding  to  Mode  I  and  Mode  II  loading  are  energy 
orthogonal,  we  have: 


Y\{ui  -l-  W//)  =  n(r//) -l- n(M//) , 

where  ui  and  Uu  are  solutions  of  the  Mode  I  and  Mode  II  loadings,  respectively.  Therefore  in 
the  case  of  combined  loading  we  have: 


G  =  Gj  +  Gjj 


4E 


(Kf+Kl). 


Closed-form  vs.  Numerical  approximation  of  G 


It  can  be  observed  from  equations  (1)  and  (4)  that  the  stresses  at  the  crack  tip  are  infinity.  It  is 
clear  that  in  a  numerical  approximation  requiring  the  integration  of  these  stresses,  even  when  the 
numerical  error  can  be  arbitrarily  reduced  by  using  an  appropriated  mesh  refinement,  this  infinite 
value  cannot  be  achieved.  Furthermore,  the  singular  behavior  of  the  stresses  at  the  crack  tip  will 
affect  the  numerical  solution  near  the  crack  tip,  and  therefore  the  quality  of  the  computed  value 
for  G  will  deteriorate  when  Aa  approaches  zero.  Because  the  quality  of  the  solution  near  the 
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crack  tip  will  affect  the  total  strain  energy  release  rate,  one  approach  would  be  to  exclude  a 
region  near  the  crack  tip  from  the  computation  of  the  energy  release  rate  (G).  This  approach  has 
limitations  however. 

Consider  the  error  in  the  computation  of  the  integral  expression  on  equation  (3)  when  excluding  a 
small  region  of  size  e  in  the  neighborhood  of  the  crack  tip  {x|  0  <  x  <  s},  that  is 


AW  =  C 


Vo  »  £ 


(Aa  -X 


dx 


X 


-  C,Aa  — 


Then  we  can  write: 


(6) 


£  rr  ^a 

r  Aa-x  .  r 

)  )  1 


'Aa-x 


dx 


Vo 


=  c,^ 


^  ^W  C,  r  Aa-x^  C  f  Aa-x^  ^  n 

G  = - =  —  J - dx  +  —  J - dx  =  C,— 

Aa  Aa  J  V  V  J  V  '  - 


X 


Aa 


X 


(7) 


Now  if  we  compute  the  second  integral  in  (7)  as  our  approximation  for  G,  we  have 


/Aa-x 

A/7  j 


Aa  V  X 


dx 


(8) 


Where  the  error  will  be  given  by 


Cj  r  Aa-x 


e  s — — —dx 
Aai\  X 


(9) 


And  the  relative  error  will  be 


Cj  r  Aa-x 


dx 


Aa-x 


C,f. 


TtAa 


dx  100% 


Considering  the  ratio  elAa  in  the  interval  [0,1]  we  can  compute  the  relative  error  associated  with 
the  first  integral  term  in  equation  (7)  as  shown  in  Figure  16. 
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Figure  16:  Relative  error  associated  with  the  use  of  equation  (8)  in  the  interval  [0, 1]  (left).  Detail  in 

the  interval  [0,  0.01]  (right). 


It  can  be  observed  that  the  contribution  of  a  small  region  near  the  crack  tip  in  the  computation  of 
the  energy  release  rate  (G)  is  of  considerable  significance.  For  example,  excluding  a  region  of  the 
crack  tip  of  0.2%  the  size  of  the  crack  increment  (s^Aa=0.002),  the  value  of  G  will  be  off  by  6%.  It 
follows  that  in  order  to  achieve  a  good  approximation  of  G,  the  stresses  near  the  crack  tip  have  to 
be  known  with  high  accuracy,  which  is  not  possible  to  achieve  with  the  finite  element  method 
without  the  use  of  highly  refined,  (geometrically  graded)  meshes  leading  to  a  significant  increase 
in  computational  cost. 


An  alternative  approach  is  discussed  in  the  following.  The  exact  solution  near  the  crack  tip  is  of 
the  form  <T  ~  Ar^~^ ,  where  /A  is  a  constant  and  A  is  the  first  eigenvalue  corresponding  to  the 
solution  of  the  corresponding  elasticity  problem.  Therefore,  near  the  crack  tip  and  can  be 
approximated  as: 


=  Ar^-' 


Then,  computing  or  T^from  the  finite  element  solution  at  two  points  (Figure  17),  it  is  possible 
to  find  A,  B,  and  A  as  follows: 
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<71 

<72 


Therefore  we  can  estimate  the  value  of  e  in  equation  (9)  by  using: 


(10) 


jcr;'’V/”’dx 


Where  w"*"”  is  obtained  from  the  numerical  approximation,  and  (7^^^  is  computed  by  evaluating 
equation  (10)  at  the  corresponding  points.  Since  equation  (11)  has  to  be  computed  numerically, 
we  cannot  include  zero  as  the  lower  limit,  therefore  we  use  the  parameter which  indicates  the 
first  point  in  the  integration  interval  neglecting  the  origin.  It  follows  that  the  quality  of  the  correction 
will  depend  on  the  number  of  points  used  for  the  extraction  of  u"™  . 

Finally  we  compute  a  corrected  value  for  G  as: 


Gc  =  G  +  e 


This  methodology  will  be  illustrated  with  examples  in  the  following. 
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Example  1  (Mode  I  loading):  Double  cantilever  beam 

The  double  cantilever  beam  shown  in  Figure  18  is  made  of  an  isotropic  material  and  loaded  in 
pure  Mode  I.  In  this  case  G=J,  where  J  is  the  J-integral.  Therefore  we  will  consider  the  computed 
value  of  J  as  the  reference  value^  for  G.  Material  properties:  E  =  1 .010e+007  psi,  v  =  0.365  (plane 
stress). 


Figure  18:  Double  cantilever  beam. 


Mesh 

The  mesh  used  for  computing  the  numerical  approximations  was  designed  so  as  to  minimize  the 
pollution  error  caused  by  the  singularity  at  the  crack  tip.  Five  layers  of  geometrically  graded 
elements  towards  the  crack  tip  where  used  as  shown  in  Figure  19.  Because  of  symmetry,  only 
one  half  of  the  beam  was  considered  for  the  analysis. 


Figure  19:  Mesh  and  boundary  conditions  for  example  1. 


Extraction 

For  each  interval  Aa  the  approximated  value  of  the  energy  release  rate  ( G  )  was  computed 
numerically  using  expression  (12),  and  analytically,  using  expression  (8). 


^  The  J-integral  is  computed  using  the  finite  element  software  StressCheck. 
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Aa 

Gi=^ - - -  (12) 

Aa 


Reference  value  (J):  The  finite  element  solution  was  obtained  by  p-extension  using  the  mesh 
shown  in  Figure  19.  The  polynomial  order  of  the  elements  was  increased  from  1  to  8  and  the 
fracture  mechanics  parameters  (Kl,  Kll  and  J)  were  extracted  for  each  one  of  the  eight  solutions 
and  the  results  shown  in  Table  3  as  a  function  of  the  number  of  degrees  of  freedom  (DOF).  It  can 
be  seen  that  the  value  of  J  remains  practically  unchanged  for  the  last  three  solutions.  The  value 
J=0.034592  was  used  as  the  reference  solution  for  G. 


Table  3:  Computed  values  of  K/,  Kn  and  J. 


DOF 

Kl 

Kll 

J 

68 

507.24 

0 

1 .0665e-02 

192 

646.05 

0 

2.8653e-02 

324 

586.50 

0 

3.1507e-02 

512 

585.97 

0 

3.3492e-02 

756 

591.85 

0 

3.4389e-02 

1056 

591.23 

0 

3.4585e-02 

1412 

590.63 

0 

3.4606e-02 

1824 

590.80 

0 

3.4592e-02 

Energy  release  rate  (G):  The  computation  of  the  approximated  energy  release  rate  G  was 
performed  using  equation  (12)  for  three  values  of  Aa,  with  the  values  of  the  stress  (Oy)  and 
displacement  {Uy)  obtained  from  the  finite  element  solution  corresponding  to  1824  DOF.  The 

J  —  G 

results  are  shown  in  Table  4.  The  percent  relative  error  was  computed  as:  - - x  100 


Table  4:  Computed  values  of  G  and  the  relative  error  (%)  with  respect  to  J. 


Aa 

e 

G 

er% 

8.1818e-03 

3.0375e-04 

2.5994e-02 

24.8  % 

1 .6363e-02 

3.0375e-04 

2.8383e-02 

17.9% 

2.4545e-02 

3.0375e-04 

2.9429e-02 

14.9% 

The  computation  of  the  correction  for  the  approximated  energy  release  rate  was  performed  by 
extracting  u""'”  a\  100  (equidistant)  points,  and  500  (equidistant)  points  in  the  interval  [0,  £].  The 

effect  of  the  addition  of  the  correction  term  in  the  value  of  G  is  shown  in  Figure  20  as  a  function 
of  Aa.  The  reference  solution  is  independent  of  Aa,  and  the  other  three  solutions  shown  in  the 
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figure  are  for  the  cases  of  no  correction  (Go),  and  with  corrections  for  100  (Gioo)  and  500  (G500) 
extraction  points.  The  relative  error  for  each  case  is  shown  in  Figure  21.  The  number  of  points  to 

extract  determines  the  value  for  that  is  going  to  limit  the  quality  of  the  approximation.  The 

results  obtained  for  100  and  500  points  are  shown  in  Tables  5  and  6. 

Table  5:  Corrected  values  of  G  and  its  relative  error  with  respect  to  Jfor  100  points. 


Aa 

^0 

Gc 

e,% 

8.1818e-03 

3.01  e-06 

3.3738e-02 

2.47  % 

1 .6363e-02 

3.01  e-06 

3.3983e-02 

1 .76  % 

2.4545e-02 

3.01  e-06 

3.4094e-02 

1 .44  % 

Table  6:  Corrected  values  of  G  and  its  relative  error  with  respect  to  Jfor  500  points. 


Aa 

^0 

Gc 

e,% 

8.1818e-03 

6.1e-07 

3.4206e-02 

1.11  % 

1 .6363e-02 

6.1e-07 

3.4323e-02 

0.  78  % 

2.4545e-02 

6.1e-07 

3.4378e-02 

0.  62  % 

Figure  20:  Computed  value  of  G  (Go)  and  its  corrections  for  100  (Gwo)  and  500  (G500)  points 


extraction. 
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Figure  21:  Computed  numerical  and  analytical  error  for  G  (Go)  and  its  corrections  for  100  (G^o)  and 

500  (Gsoo)  points  extraction. 


Correction  using  and  extrapolation  function  for  o  and  u:  The  results  shown  in  Tables  5  and  6 
indicate  that  the  solution  still  depends  on  the  value  of  since  the  numerical  integral  for  the 
correction  term  must  exclude  the  origin.  An  alternative  approach  was  considered  to  remove  this 
limitation.  The  problem  in  using  equation  (11)  is  that  the  integral  has  to  be  performed  numerically 
because  the  displacement  is  known  from  the  numerical  solution,  and  that  precludes  extending  the 
integration  limit  to  include  zero.  The  solution  is  to  perform  the  integral  analytically  which  requires 
an  approximation  of  the  displacement  function  in  the  integration  interval.  Near  the  crack  tip  the 

displacement  Uy  is  of  the  form  =  Br^  .  Therefore  it  is  possible  to  write  an  approximation  for  Uy 


as: 


3^1 


where  1  and  2  refer  to  two  points  along  the  crack  tip  coordinate  system  as  shown  in  Figure  17. 
The  correction  for  the  energy  density  function  (G)  using^: 


^  The  reason  for  using  and  X2,  instead  of  a  single  value  X  is  because  these  values  are  computed  from 
fitting  numerical  data,  and  therefore  they  are  expected  to  be  slightly  different,  whereas  in  theory  they  should 
be  the  same. 
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£  £ 

J  Jx  AB^  x^‘"‘  (Aa  -  x)^  dx 


Aa  Aa 


With  A  given  in  equation  (10).  This  integral  can  be  written  as: 

A  6  (1  +  A).  k\ 

Where  any  term  between  brackets,  such  as  (/Ij  ,  is  expanded  as  (/Ij  -1)  ■ 

We  refer  to  Mathematica.'^  We  observe  that  the  solution  for  this  integral  exist  and  is  bounded, 
however  we  need  to  compute  a  series  expansion  in  order  to  obtain  an  approximated  value. 
Alternatively,  since  the  function  Uy  is  not  extracted  at  the  crack  tip  but  at  the  end  of  the  crack 
increment  Aa  were  the  function  u  is  very  smooth,  we  can  approximate  Uy  by  a  linear  function,  that 
is  Uy  =  Bx+C.  Then, 


It  follows  that 


B  = 


yi 


C  =  -Bx, 

yi  ^ 


£  £ 

I  fT^g^r/x  I  Ax^-'  [B{Aa  -  x)  +  C\lx 


Aa  Aa 


Aa 


t  t 

■  B^  x^dx  +  (C  +  5Ag)|  x^~^dx 


Aa 


„/l+l  „/l 

-B - +  {C  +  BAa) — 

/I  + 1  X 


(13) 


The  case  for  mode  II  loading  is  analogous.  The  results  shown  in  Table  7  were  obtained  using 
equation  (13)  to  compute  the  correction  forG  . 


Mathematica  is  a  trademark  of  Wolfram  Research,  Inc. 
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Table  7:  Values  of  Gc  and  its  relative  error  with  respect  to  J  using  an  approximated  closed  form  for 

the  correction  term. 


Aa 

J 

Gc 

er% 

8.1818e-03 

3.4592e-02 

3.4596e-02 

0.0125% 

1.6363e-02 

3.4592e-02 

3.4605e-02 

0.0388  % 

2.4545e-02 

3.4592e-02 

3.461 6e-02 

0.0681  % 

Example  2  (Mode  II  loading):  Plate  with  a  crack  under  pure  shear. 


S 

SiF 


s 


s 

Figure  24:  Plate  with  a  crack  under  pure  shear. 


The  problem  shown  in  Figure  24  corresponds  to  a  flat  panel  of  constant  thickness  with  a  central 
crack  under  pure  shear.  Isotropic  material  properties  are  considered  in  order  to  provide  a 
reference  solution.  Since  this  case  is  in  pure  Mode  II  the  value  of  the  energy  release  rate  (G/;)  will 
be  coincident  with  the  value  of  the  J  integral  {J).  Therefore  we  will  consider  the  computed  value  of 
J  as  the  reference  value  for  G/;.  Material  properties:  £  =  1 .010e+007  psi,  v  =  0.365  (plane  stress). 


Mesh 


The  mesh  used  for  computing  the  numerical  approximations  was  designed  so  as  to  minimize  the 
pollution  error  caused  by  the  singularity  at  the  crack  tip.  Four  geometrically  graded  layers  of 
elements  where  used  with  such  purpose  as  shown  in  Figure  25.  Because  of  symmetry,  one 
quarter  of  the  plate  was  considered  for  the  analysis. 
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Extraction 

For  each  interval  Aa  the  approximated  value  (G  )  of  the  energy  release  rate  (G)  was  computed 
both  numerically,  using  expression  (14),  and  analytically,  using  expression  (8). 

Aa 

Gu=^ - - -  (14) 

Aa 


Reference  value  (J):  The  reference  solution  was  obtained  by  p-extension  using  the  finite 
element  mesh  shown  in  Figure  25.  The  polynomial  order  of  the  elements  was  increased  from  1  to 
8  and  the  fracture  mechanics  parameters  were  extracted  for  each  one  of  the  eight  solutions  and 
the  results  shown  in  Table  8  as  a  function  of  the  number  of  degrees  of  freedom  (DOF).  It  can  be 
seen  that  the  value  of  J  remains  practically  unchanged  for  the  last  three  solutions.  The  value 
J=0. 39661  was  used  as  the  reference  solution. 


Table  8:  Computed  values  of  Ki,  Kn  and  J. 


DOF 

KI 

Kll 

J 

59 

0 

-1428.7 

1.8779e-01 

171 

0 

-1913.3 

3.3865e-01 

291 

0 

-2031.0 

3.9681  e-01 

463 

0 

-2014.6 

4.0095e-01 

687 

0 

-2001.5 

3.991 1e-01 

963 

0 

-1998.2 

3.9723e-01 

1291 

0 

-1999.6 

3.9655e-01 

1671 

0 

-2001.5 

3.9661e-01 
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Energy  release  rate  (G):  The  computation  of  the  approximated  energy  release  rate  G  was  done 
using  equation  (14)  for  three  values  of  Aa,  with  the  values  of  the  stress  (Xxy)  and  displacement 
{Ux)  obtained  from  the  finite  element  solution  corresponding  to  1671  DOF.  The  results  are  shown 
in  Table  9. 


Table  9:  Computed  values  of  G  and  its  relative  error  with  respect  to  J. 


Aa 

e 

G 

er% 

2.25e-03 

5.0625e-04 

1.6645e-01 

58.12% 

4.50e-03 

5.0625e-04 

2.3068e-01 

41.96% 

6.75e-03 

5.0625e-04 

2.6037e-01 

34.49  % 

The  computation  of  the  correction  for  the  approximated  energy  release  rate  was  performed  using 
expression  analogous  to  equation  (13),  and  the  results  are  shown  in  Table  10. 


Table  10:  Corrected  values  of  G  and  its  reiative  error  with  respect  to  J  using  an  approximated  ciosed 


form  for  e. 


Aa 

J 

Gc 

er%  (comd) 

2.25e-03 

3.9661e-01 

3.9906e-01 

0.402  % 

4.50e-03 

3.9661e-01 

3.9704e-01 

0.107% 

6.75e-03 

3.9661e-01 

3.9682e-01 

0.161  % 

Example  3  (Mode  I  loading) 

Consider  the  same  model  problem  of  example  1  (Figure  18)  using  orthotropic  materials  properties 
(Ell  =  8.0e6  psi,  Ejt  =  I.OeG  psi,  Vlt  =  0.3,  Glt  =  0.6e6  psi,  plane  stress).  Since  there  is  material, 
geometric  and  load  symmetry  the  problem  remains  a  pure  Mode  I  loading  and  mesh  shown  in 
Figure  19  was  used  for  the  analysis.  Therefore  we  can  use  J  as  the  reference  value  for  G;. 


Reference  value  (J):  The  reference  solution  was  obtained  by  p-extension  and  the  results  shown 
in  Table  1 1 .  The  value  J=0.052817  was  used  as  the  reference  solution. 

Table  11:  Computed  values  ofJ. 


DOF 

J 

68 

2.7688e-02 

192 

4.601 9e-02 

324 

4.9508e-02 

512 

5.2128e-02 

756 

5.2759e-02 

1056 

5.2845e-02 

1412 

5.2843e-02 

1824 

5.281 7e-02 
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Energy  release  rate  (G):  The  computation  of  the  approximated  energy  release  rate  G  was  done 
using  equation  (12)  with  the  values  of  the  stress  (oy)  and  displacement  (t/y)  obtained  from  the 
finite  element  solution  corresponding  to  1824  DOF.  The  results  are  shown  in  Table  12. 

Table  12:  Computed  values  of  G  and  its  relative  error  with  respect  to  J. 


Aa 

£ 

G 

er% 

8.1819e-03 

3.038e-04 

3.9780e-02 

24.65  % 

1.6364e-02 

3.038e-04 

4.3478e-02 

17.68% 

2.4545e-02 

3.038e-04 

4.5102e-02 

14.61  % 

Using  equation  (13)  to  compute  the  correction  for  G,  the  results  shown  in  Table  13  were 
obtained,  providing  a  clear  indication  that  the  benefits  of  the  approach  in  the  computation  of  the 
energy  release  rate  demonstrated  for  isotropic  materials  is  also  realized  for  orthotropic  materials. 


Table  13:  Values  of  G  and  its  reiative  error  with  respect  to  J  using  an  approximated  ciosed  form  for 

the  correction  term. 


Aa 

J 

Gc 

er% 

8.1819e-03 

5.281 7e-02 

5.281 7e-02 

1.41e-04% 

1 .6364e-02 

5.281 7e-02 

5.281 6e-02 

1 .08e-03  % 

2.4545e-02 

5.281 7e-02 

5.2821  e-02 

6.61e-03% 

Example  4  (Mode  II  loading) 

Consider  once  again  the  model  problem  of  example  2  (Figure  24)  however  this  time  using 
orthotropic  materials  properties 

(Ell  =  8.0e6  psi,  Ejt  =  1.0e6  psi,  Vlt  =  0.3,  Glt  =  0.6e6  psi,  plane  stress).  Since  there  is  material, 
geometric  symmetry  and  load  anti-symmetry,  the  plate  is  under  pure  Mode  II  loading  when  using 
an  orthotropic  material  aligned  with  the  geometric  axis.  Therefore  we  can  use  J  as  the  reference 
value  for  G/;. 

Reference  value  (J):  The  reference  solution  was  obtained  by  p-extension  and  the  results  shown 
in  Table  14.  The  value  J=1. 0315  was  used  as  the  reference  solution. 
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Table  14:  Computed  values  of  J. 


DOF 

J 

68 

2.5843e-02 

197 

6.5063e-01 

334 

8.2705e-01 

531 

9.3559e-01 

788 

9.9751  e-01 

1105 

1.0204e+00 

1482 

1.0287e+00 

1919 

1.0315e+00 

Energy  release  rate  (G):  Computation  of  the  approximated  energy  release  rate  G  was  done  as 
described  before  and  the  results  are  shown  in  Table  15. 

Table  15:  Computed  values  of  G  and  its  relative  error  with  respect  to  J. 


Aa 

£ 

G 

er% 

2.25e-03 

5.0625e-04 

4.3204e-01 

58.12% 

4.50e-03 

5.0625e-04 

5.9989e-01 

41.84  % 

6.75e-03 

5.0625e-04 

6.7736e-01 

34.33  % 

Using  equation  (12)  to  compute  the  correction  for  Gin,  we  obtained  the  results  shown  in  Table 
16.  Note  that  for  this  case  the  error  was  still  under  1%  but  larger  than  for  the  other  model 
problems. 


Table  16:  Values  of  G  and  its  relative  error  with  respect  to  J  using  an  approximated  closed  form  for 


the  correction  term. 


Aa 

J 

Gc 

er% 

2.25e-03 

1.0315e+00 

1.0254e+00 

0.59  % 

4.50e-03 

1.0315e+00 

1.0247e+00 

0.66  % 

6.75e-03 

1.0315e+00 

1.0261e+00 

0.53  % 

Example  5  (Mixed  mode):  Plate  with  a  crack  under  shear  and  axial 
tension. 

In  order  to  provide  a  case  of  study  for  a  mixed  mode  loading,  and  at  the  same  time  have  enough 
information  to  correlate  results  with  other  solution  techniques,  the  problem  shown  in  Figure  27 
was  considered,  representing  a  flat  panel  of  constant  thickness  with  a  central  crack  subjected  to 
the  combined  effects  of  shear  and  tension  loads.  The  material  properties  are  the  same  as  those 
used  for  problems  3  and  4. 
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Figure  27:  Proposed  model  for  mixed  mode. 


Since  in  this  example  Gu  has  to  be  equal  than  the  one  computed  for  Example  4,  and  Gt=  Gi+  Gu 
=  J,  we  can  use  this  values  as  a  reference  to  determine  the  accuracy  of  the  method  under  a 
mixed  mode. 


Mesh 

The  mesh  used  for  computing  the  numerical  approximations  was  designed  so  as  to  minimize  the 
pollution  error  caused  by  the  singularity  at  the  crack  tip.  Four  geometrically  graded  layers  of 
elements  where  used  with  such  purpose  as  shown  in  Figure  28. 


Figure  28:  Mesh  and  boundary  conditions  for  exampie  5. 
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Reference  value  (J):  The  reference  solution  was  obtained  by  p-extension  using  the  finite 
element  mesh  shown  in  Figure  28.  The  polynomial  order  of  the  elements  was  increased  from  1  to 
8  and  the  J-integral  was  extracted  for  each  one  of  the  eight  solutions  and  the  results  shown  in 
Table  17  as  a  function  of  the  number  of  degrees  of  freedom  (DOF).  It  can  be  seen  that  the  value 
of  J  remains  practically  unchanged  for  the  last  three  solutions.  The  value  J=3.9865  was  used  as 
the  reference  solution. 


Table  17:  Computed  values  of  J. 


DOF 

J 

255 

2.6171e+00 

737 

3.1841e+00 

1251 

3.5666e+00 

1989 

3.7891  e+00 

2951 

3.9091  e+00 

4137 

3.9606e+00 

5547 

3.9800e+00 

7181 

3.9865e+00 

Reference  value  (G//):  Assuming  superposition,  the  reference  solution  for  the  energy  release 
rate  corresponding  to  the  second  (shearing)  mode  Gn  is  provided  by  the  energy  release  rate 
computed  in  Example  4,  while  the  value  of  Jean  be  compared  with  the  sum  of  G/  +  G//.  Table  18 

shows  the  results  of  the  computation  of  the  approximated  energy  release  rate  G  utilizing 
equations  (12)  and  (14). 

Table  18:  Computed  values  of  Gi ,  Gn ,  Gt  =  Gi  Gn  and  its  relative  error  with  respect  to  J. 


Aa 

e 

Gi 

Gn 

Gt 

J 

er% 

2.25e-03 

5.0625e-04 

1.2378 

4.3241e-01 

1.6702 

3.9865 

58.10% 

41.89% 

34.40  % 

4.50e-03 

5.0625e-04 

1 .7246 

5.9189e-01 

2.3165 

3.9865 

6.75e-03 

5.0625e-04 

1 .9526 

6.6254e-01 

2.6151 

3.9865 

Using  equation  (13)  to  compute  the  correction  for  Gi  andGn ,  the  results  shown  in  Tables  19 
and  20  were  obtained. 


Tabie  19:  Corrected  vaiues  of  Gi ,  Gn ,  Gt  and  its  reiative  error  with  respect  to  J  using  an 
approximated  ciosed  form  for  the  correction  term. 


Aa 

GfC 

G//C 

GfC 

J 

er% 

2.25e-03 

2.9561 

1.0157 

3.9718 

3.9865 

0.37  % 

4.50e-03 

2.9650 

1.0033 

3.9683 

3.9865 

0.46  % 

6.75e-03 

2.9739 

9.97e-01 

3.9710 

3.9865 

0.39  % 
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Table  20:  Corrected  values  of  Gn  and  its  relative  error  with  respect  to  Jref  obtained  in  Example  4. 


Aa 

G//C 

Jref 

er% 

2.25e-03 

1.0157 

1.0315 

1 .54  % 

4.50e-03 

1.0033 

1.0315 

2.73  % 

6.75e-03 

9.97e-01 

1.0315 

3.45  % 

It  is  observed  that  the  computation  of  G/  and  Gn  for  a  mixed  mode  loading  provides  results  with 
acceptable  error  bounds.  The  method  is  considered  to  be  adequate  for  the  computation  of  G/  and 
G//  for  composite  materials  in  2D.  The  extension  of  the  proposed  methodology  for  the 
computation  of  the  energy  release  rate  components  in  three-dimensions  will  be  included  in  the 
proposal  for  the  Phase  II  project. 

Summary  and  Conclusions 

One  of  the  main  goals  of  this  project  was  the  investigation  of  mathematical  models  for  the 
simulation  of  structural  stability  and  post-buckling  responses  of  light-weight  unitized  airframe 
components  fabricated  from  aluminum  plates  by  high  speed  machining  techniques.  It  has  been 
observed  in  machining  experiments  that  some  thin-walled  components  exhibit  local  buckling  even 
in  the  unloaded  condition,  due  to  residual  stresses.  Therefore  the  effects  of  residual  stresses 
were  taken  into  account  in  simulations  performed  for  the  purposes  of  computing  the  buckling 
behavior  with  and  without  external  loads. 

The  effect  of  residual  stresses  on  the  structural  stability  of  thin  unitized  components  machined 
from  aluminum  plates,  in  particular  7050-T74  and  7050-T7451  plates,  was  investigated.  The 
findings  indicated  that  residual  stresses  induced  in  the  plate  during  the  rolling  operation  (bulk 
stresses)  and  residual  stresses  caused  by  high  speed  machining  (machining-induced  residual 
stresses)  have  a  significant  effect  on  the  stability  and  should  be  included  as  part  of  the  modeling 
assumption  when  designing  thin  unitized  components.  It  was  also  found  that  bulk  residual 
stresses  are  not  large  enough  to  cause  the  specimen  to  buckle  after  removing  the  machined  part 
from  the  fixture,  but  they  may  have  a  detrimental  effect  on  the  buckling  load  once  the  component 
is  loaded.  Machining-induced  residual  stresses  on  the  other  hand  are  capable  of  producing 
buckling  once  the  part  is  removed  from  the  fixture. 

The  second  goal  was  to  develop  a  procedure  to  determine  the  residual  strength  of  structures 
made  of  composite  materials  damaged  by  delamination.  To  that  goal,  the  implementation  of  a 
modified  virtual  crack  closure  technique  for  the  computation  of  the  energy  release  rate  in  two 
dimensions  was  investigated.  The  advantages  of  the  algorithm  was  demonstrated  by  solving 
model  problems  in  Mode  I,  Mode  II  and  mixed  mode  loading  for  isotropic  and  orthotropic 
materials. 
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Future  work 

Residual  stresses  and  structural  stability 

The  preliminary  procedure  developed  to  incorporate  the  effects  of  residual  stresses,  either 
material  (bulk)  or  machining-induced,  in  the  onset  of  stability  of  unitized  components  was 
demonstrated  to  be  highly  effective.  However,  the  procedure  is  very  time-consuming  and  needs 
to  be  optimized  and  automated.  There  are  also  some  limitations  on  the  current  implementation 
that  need  to  be  removed.  Since  the  residual  stresses  are  introduced  by  means  of  a  thermal  load, 
the  application  of  external  loads  is  restricted  to  loads  that  satisfy  the  equilibrium  condition,  since 
additional  constraints  can  affect  the  desired  residual  stress  distribution.  These  activities  will  be 
incorporated  as  part  of  the  Phase  II  proposal. 

Delamination  of  composite  materials 

The  modified  virtual  crack  closure  technique  investigated  during  the  Phase  I  project  provides  an 
effective  and  reliable  procedure  for  the  computation  of  the  strain  energy  release  rates  for  mixed 
modes  in  orthotropic  materials.  Research  that  will  address  implementation  of  the  methodology 
within  the  framework  of  StressCheck  and  its  extension  to  three-dimensional  analysis  will  be  one 
of  the  main  objectives  of  the  Phase  II  project. 
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